#install.packages("ISwR")
library(ISwR)

# 4. Descriptive statistics and graphics ----------------------------------


# 4.1 Summary statistics for a single group -------------------------------

x <- rnorm(50)
mean(x)
sd(x)
var(x)
median(x)
quantile(x)
quantile(x,seq(0,1,0.1))

attach(juul)
mean(igf1) # si hay NAs mean no funciona
mean(igf1, na.rm=T)  # ahora si

length(igf1) # length tiene en cuenta los NAs
sum(!is.na(igf1)) # asi contamos solo los que no son NAs

summary(igf1)
summary(juul)

# creamos los factores
detach(juul)
juul$sex <- factor(juul$sex,labels=c("M","F"))
juul$menarche <- factor(juul$menarche,labels=c("No","Yes"))
juul$tanner <- factor(juul$tanner, labels=c("I","II","III","IV","V"))
attach(juul)
summary(juul)

# o de otra manera mas compacta
juul <- transform(juul,
                  sex=factor(sex,labels=c("M","F")),
                  menarche=factor(menarche,labels=c("No","Yes")),
                  tanner=factor(tanner,labels=c("I","II","III","IV","V")))


# 4.2 Graphical display of distributions ----------------------------------


# 4.2.1 Histograms --------------------------------------------------------

hist(x)

# la x del histograma no tiene porque ser constante
mid.age <- c(2.5,7.5,13,16.5,17.5,19,22.5,44.5,70.5)
acc.count <- c(28,46,58,20,31,64,149,316,103)
age.acc <- rep(mid.age,acc.count)

hist(age.acc,breaks=c(0,5,10,16,17,18,20,25,60,80))
# the area of a column is proportional to the number. The y-axis is in density units
# (that is, proportion of data per x unit), so that the total area of the histogram will be 1.

# 4.2.2 Empirical cumulative distribution ---------------------------------

n <- length(x)
plot(sort(x),(1:n)/n,type="s",ylim=c(0,1))

# 4.2.3 Q–Q plots ---------------------------------------------------------

qqnorm(x)

# 4.2.4 Boxplots ----------------------------------------------------------

par(mfrow=c(1,2))
boxplot(IgM)
boxplot(log(IgM))
par(mfrow=c(1,1))

# 4.3 Summary statistics by groups ----------------------------------------

attach(red.cell.folate)
red.cell.folate
tapply(folate,ventilation,mean)
tapply(folate,ventilation,sd)
tapply(folate,ventilation,length)
# for a nicer display
xbar <- tapply(folate, ventilation, mean)
s <- tapply(folate, ventilation, sd)
n <- tapply(folate, ventilation, length)
cbind(mean=xbar, std.dev=s, n=n)
#podemos pasarle parametros a la funcion
tapply(igf1, tanner, mean, na.rm=T)

aggregate(juul[c("age","igf1")], list(juul$tanner,sex=juul$sex), mean, na.rm=T)

by(juul, juul["sex"], summary)

# 4.4 Graphics for grouped data -------------------------------------------


# 4.4.1 Histograms --------------------------------------------------------

attach(energy)
expend.lean <- expend[stature=="lean"]
expend.obese <- expend[stature=="obese"]
par(mfrow=c(2,1))
hist(expend.lean,breaks=10,xlim=c(5,13),ylim=c(0,4),col="white")
hist(expend.obese,breaks=10,xlim=c(5,13),ylim=c(0,4),col="grey")
par(mfrow=c(1,1))

# 4.4.2 Parallel boxplots -------------------------------------------------

boxplot(expend ~ stature) # y ~ x should be read “y described using x”
boxplot(expend.lean,expend.obese)

# 4.4.3 Stripcharts -------------------------------------------------------

opar <- par(mfrow=c(2,2), mex=0.8, mar=c(3,3,2,1)+.1)
stripchart(expend ~ stature)
stripchart(expend ~ stature, method="stack")
stripchart(expend ~ stature, method="jitter")
stripchart(expend ~ stature, method="jitter", jitter=.03)
par(opar)

# 4.5 Tables --------------------------------------------------------------


# 4.5.1 Generating tables -------------------------------------------------

# creamos una tabla
caff.marital <- matrix(c(652,1537,598,242,36,46,38,21,218,327,106,67),nrow=3,byrow=T)
caff.marital
colnames(caff.marital) <- c("0","1-150","151-300",">300")
rownames(caff.marital) <- c("Married","Prev.married","Single")
caff.marital
names(dimnames(caff.marital)) <- c("marital","consumption")
caff.marital
as.data.frame(caff.marital)

table(sex)
table(menarche,tanner)
table(sex,menarche)

table(tanner,sex)
xtabs(~ tanner + sex, data=juul)

#cuando tenemos mas dimensiones:
ftable(coma + diab ~ dgn, data=stroke)

# 4.5.2 Marginal tables and relative frequency ----------------------------

tanner.sex <- table(tanner,sex)
tanner.sex
margin.table(tanner.sex,1)
margin.table(tanner.sex,2)
table(tanner)
table(sex)

prop.table(tanner.sex,1) # para mostrar las frecuencias relativas
round(prop.table(tanner.sex,1)*100, 1) # en %

round((tanner.sex/sum(tanner.sex))*100,1) # %s del total

# 4.6 Graphical display of tables -----------------------------------------


# 4.6.1 Barplots ----------------------------------------------------------

total.caff <- margin.table(caff.marital,2)
total.caff
barplot(total.caff, col="white")

par(mfrow=c(2,2))
barplot(caff.marital, col="white")
barplot(t(caff.marital), col="white")
barplot(t(caff.marital), col="white", beside=T)
barplot(prop.table(t(caff.marital),2), col="white", beside=T)
par(mfrow=c(1,1))

barplot(prop.table(t(caff.marital),2),beside=T, 
        legend.text=colnames(caff.marital), 
        col=c("white","grey80","grey50","black"))


# 4.6.2 Dotcharts ---------------------------------------------------------

# contains the same information as barplots with beside=T but give quite a different visual impression
dotchart(t(caff.marital), lcolor="black")

# 4.6.3 Piecharts ---------------------------------------------------------

opar <- par(mfrow=c(2,2),mex=0.8, mar=c(1,1,2,1))
slices <- c("white","grey80","grey50","black")
pie(caff.marital["Married",], main="Married", col=slices)
pie(caff.marital["Prev.married",],main="Previously married", col=slices)
pie(caff.marital["Single",], main="Single", col=slices)
par(opar)



# 4.7 Exercises -----------------------------------------------------------

# 4.1 Explore the possibilities for different kinds of line and point plots.
# Vary the plot symbol, line type, line width, and colour.
plot(sort(rnorm(25)), col="red", pch=3, type="b", lwd=3)

# 4.2 If you make a plot like plot(rnorm(10),type="o") with overplotted
# lines and points, the lines will be visible inside the plotting
# symbols. How can this be avoided?
plot(rnorm(10),type="o", pch=21, bg="white" )
plot(1:25, pch=1:25)
text((1:25)+1, 1:25)

# 4.3 How can you overlay two qqnorm plots in the same plotting area?
# What goes wrong if you try to generate the plot using type="l", and
# how do you avoid that?
data1 <- rnorm(100)
data2 <- rnorm(100)
res1 <- qqnorm(sort(data1), plot.it=F)
res2 <- qqnorm(sort(data2), plot.it=F)
rangex <- range(res1$x, res2$x)
rangey <- range(res1$y, res2$y)
plot(res1, xlim=rangex, ylim=rangey, type="n")
points(res1$x,res1$y, col="red", type="l")
points(res2, col="green", type="l")

# 4.4 Plot a histogram for the react data set. Since these data are highly
# discretized, the histogram will be biased. Why? You may want to try
# truehist from the MASS package as a replacement.
react
summary(react)
table(react)
library(MASS)

opar <- par(mfrow=c(1,2))
hist(react, breaks=20)
truehist(react,h=1,x0=.5)
par(opar)
# 4.5 Generate a sample vector z of five random numbers from the uniform
# distribution, and plot quantile(z,x) as a function of x (use curve, for instance).
z <- runif(5)
data <-sapply(seq(5, 205, by=50), function (i){quantile(runif(i))})
plot(1,1, xlim=c(0,1), ylim=c(0,1),type="n")
sapply(1:5, function(i){lines(seq(0,1,by=0.25),data[,i], col=i)})







